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1  PREFACE 

This  study  was  conducted  for  the  Explosion  Effects  Division  (EED),  Structures  Labo¬ 
ratory  (SL),  of  the  U.  S.  Army  Engineer  Waterways  Experiment  Station  (WES).  The 
work  was  sponsored  by  Field  Command,  Defense  Nuclear  Agency,  under  a  Scientific 
Services  Agreement  issued  by  Battelle,  Contract  No.  DAAL03-86-D-0001,  Delivery 
Order  1308. 

The  research  was  accomplished  by  Dr.  C.  Wayne  Mastin,  Mississippi  State  Univer¬ 
sity,  under  the  technical  supervision  of  Mr.  Howard  G.  White,  EED,  and  the  general 
supervision  of  Mr.  L.  K.  Davis,  Chief,  EED,  and  Mr.  Bryant  Mather,  Director,  SL. 

At  the  time  of  publication  of  this  report,  Director  of  WES  was  Dr.  Robert  W. 
Whalin.  Commander  was  COL  Leonard  G.  Hassell,  EN. 
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AN  AUTOMATIC  ADAPTIVE  REZONE 
SCHEME  FOR  HULL  CODE 


3  INTRODUCTION 

The  success  or  failure  of  a  hydrodynamic  computation  often  depends  on  the  grid 
that  is  used.  This  is  especially  true  in  water  shock  problems  where  large  solution 
gradients  must  be  resolved  by  the  numerical  algorithm.  An  insufficient  number  of 
grid  points  in  the  region  near  a  shock  wave  leads  to  either  oscillations  or  smearing  in 
the  solution.  A  grid  with  uniform  spacing  in  each  coordinate  direction  could  be  used 
for  solving  problems  of  this  type,  but  that  would  be  wasteful  since  there  are  typically 
large  regions  where  the  solution  is  nearly  constant.  It  would  be  desireable  to  have  a 
grid  with  a  high  concentration  of  grid  points  in  regions  where  the  solution  gradients 
are  large,  and  very  few  points  where  the  solution  is  nearly  constant.  However,  this 
cannot  be  done  a  priori ,  since  the  shock  wave  location  is  determined  by  the  solution. 
What  is  needed  is  an  adaptive  rezoning  capability  that  regenerates  the  grid  at  each 
step  of  the  computational  procedure,  based  on  the  current  values  of  the  numerical 
solution. 

The  HULL  code  has  several  existing  rezone  capabilities  for  two-dimensional  prob¬ 
lems;  however,  none  of  methods  described  in  the  accompanying  documentation  [1] 
are  truly  general  purpose  solution  adaptive  algorithms.  This  report  describes  the 
development  of  an  adaptive  grid  algorithm  for  the  two-dimensional,  Eulerian,  finite- 
difference  method.  The  same  technique  could  be  used  in  three  dimensions,  but  there 
was  no  three-dimensional  rezoning  routine  in  the  existing  program.  There  was  a  mov¬ 
ing  grid  procedure  which  could  be  helpful  in  future  work  of  this  type.  The  adaptive 
rezoning  method  has  been  used  in  the  solution  of  several  problems  involving  under¬ 
water  explosions.  The  method  has  worked  well  on  the  these  problems,  given  the  grid 
restrictions  of  the  HULL  code;  namely,  the  grid  must  be  composed  of  horizontal  and 
vertical  lines. 


4  TECHNICAL  DISCUSSION 


The  HULL  code  solves  the  equations  of  continuum  mechanics  in  two  or  three  di¬ 
mensions.  For  an  axisymmetric  solution,  the  equations  can  be  written  in  cylindrical 
coordinates  as: 


Du  dp 
P~Dt  +  dx 


Dp 

lTt+P 


1  dxu  dv 

- 1 - 

x  dx  dy 


X  dx^^  +  dy^  +  xSee 


=  0 
=  0 


(1) 

(2) 
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DE 


Dt  x  dx 


PDt  +  dy  +  xdx^Sxy)  dyiSyy) 
-S-[X(uSIX  ~  vSxy  -  up)]  +  S-  -  vSyy  +  Vp) 


dy 


where 


~pg  (3) 
-pvg  (4) 


x,y  = 
u,v  = 

Sap  = 

P  = 
E  = 

P  = 
9  = 


radial  and  axial  coordinates 
radial  and  axial  velocity  components 
deviatoric  stress  tensor  components 
material  density 
total  specific  energy 
hydrostatic  pressure 
gravitational  body  force 


The  difference  scheme  used  to  solve  these  equations  is  based  on  a  method  of  Gentry, 
Martin,  and  Daly  [2].  That  method  is  itself  derived  from  the  Lagrangian  Particle- In- 
Cell  (PIC)  method  cf  Harlow  [3].  The  finite  difference  grid  is  composed  of  vertical 
and  horizontal  lines 


X  —  2  —  2  max 

y  ~  Vj,  J ~  0;  1>  '  ‘  '  i  Jmax 

The  spatial  intervals  are 

Ax,'  =  x,_  i,  i  —  1,  •  •  •  ,  Zmax 

Aj/j  =  pj  —  Pj- 1 ,  j  1 ,  •  •  •  ,  jma.x 

The  values  of  the  numerical  solution  are  defined  at  the  cell  centers  so  that,  for  a 
typical  solution  variable  /, 

hi  =  \{y}  +  Vj-i)) 

A  two-step  numerical  method  is  used  to  advance  the  solution  from  time  step  t  =  tn 
to  time  step  t  —  tn+1.  The  first  step  is  a  Lagrangian  calculation  which  assumes  the 
grid  moves  with  the  material.  This  step  advances  the  solution  using  the  above  system 
(1-4)  without  the  convection  terms.  The  next  step  is  the  transfer  of  the  Lagrangian 
solution  to  the  actual  computational  grid  at  time  fn+i  in  such  a  way  that  the  mass, 
momentum,  and  energy  are  conserved. 

An  adaptive  grid  must  sense  variations  in  the  solution  and  adjust  the  grid  accord¬ 
ingly.  Since  there  are  several  solution  variables,  more  than  one  adaptive  grid  could 
be  constructed.  However,  in  explosives  problems,  the  pressure  is  a  key  variable  in 
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assessing  blast  effects.  Thus,  the  rezoning  algorithm  will  be  based  on  the  variations 
in  the  pressure. 

Since  the  same  scheme  is  used  to  redistribute  points  in  each  coordinate  direction, 
only  the  redistribution  along  the  x-axis  will  be  explained  in  detail.  The  first  con¬ 
sideration  will  be  the  selection  of  a  weight  function  to  be  used  to  control  the  grid 
spacing.  Each  weight  function  will  control  spacing  in  only  one  direction  and  must  be 
a  function  of  a  single  variable.  Therefore  we  will  start  out  by  defining  a  function  P 
by  the  formula 

P(x{ )  =  max  {p(xi,yj) :  0  <  j  <  jmaI} 

In  order  that  the  grid  spacing  be  influenced  by  the  values  of  both  P  and  its  derivatives, 
the  weight  function  will  be  given  as 

w{Xi)  =  Co|P(x,-)|  +  c1|/>,(xi)l  +  c2\P"(Xi)\ 


The  coefficients  c,  are  included  to  add  flexibility  to  the  scheme.  If  co  is  the  dominant 
coefficient,  then  the  grid  spacing  will  be  smallest  where  the  pressure  is  greatest.  The 
C\  term  causes  grid  points  to  cluster  where  the  pressure  gradient  is  largest,  such 
as  near  shock  waves.  The  c2  term  would  cluster  points  where  large  changes  in  the 
pressure  gradient  occur,  such  as  near  oscillations  in  the  numerical  solution.  Nearly  all 
numerical  algorithms  work  best  when  there  is  a  smooth  change  in  grid  spacing.  Since 
w  is  defined  in  terms  of  P  and  its  derivatives,  it  may  change  drastically  from  point 
to  point  and  from  time  step  to  time  step.  In  such  cases  it  is  advisable  to  smooth  the 
function  P  before  computing  w.  A  diffusion-based  smoother,  given  by  the  following 
formula,  is  used  in  our  algorithm: 

P{xi)  =  P(x%)  +  T(P(xi+i)  +  P(xi- 1)  -  2 P{ii)  -  -a{P{xi+1)  -  P(x,-i))) 


where 


T 


a 


Ax,+1  +  Axi_j  +  2Ax, 


After  each  smoothing  sweep,  the  endpoint  values  are  reset  to  the  values  of  their 
interior  neighbors: 


P(x  o)  =  P{x  i) 


At  least  one  smoothing  step  is  recommended;  however,  it  should  be  noted  that  re¬ 
peated  smoothing  will  reduce  the  function  P,  and  hence  w,  to  a  constant  which  would 
result  in  a  uniform  grid. 
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The  weight  function  which  has  been  defined  is  for  rezoning  in  the  x-direction  and 
will  henceforth  be  denoted  as  wx.  If  x  is  replaced  by  y  and  i  is  replaced  by  j,  an 
analogous  weight  function  wy  can  be  derived  for  rezoning  in  the  y-direction.  Now  the 
maximum  weight  function  value  in  both  directions  can  be  computed  as 


=  max  | max  xvx(xi),  max  u>j,(y,)  j 


A  parameter  u>  >  1  will  be  introduced  which  will  determine  the  degree  to  which  the 
grid  adapts  to  the  weight  functions.  In  order  to  see  how  this  is  done,  consider  the 
grid  function 

Wi  =  i  +  (w  — (5) 
We  are  going  to  rezone  so  that  the  intervals  Ax,  for  the  new  grid  will  satisfy 

fixiWi  =  Cx  (6) 


where  Cx  is  a  constant  to  be  determined.  Now  it  can  be  seen  that  if  wx  is  near  the 
maximum  value  tCmax,  then  Wx  is  approximately  equal  to  a /  and  Ax,  ~  Cx/u>.  On  the 
other  hand,  if  wx  is  much  less  than  wmAX,  then  W,  is  nearly  1  and  Ax,  ~  Cx.  Therefore, 
when  wx  assumes  values  near  both  extremes,  the  quantity  u  is  approximately  the  ratio 
of  the  maximum  to  minimum  grid  spacing.  To  determine  the  constant  Cx ,  note  that 

*m»x  p  *max  1 

*«">«  -  X0=  Al.  =  X)  Ty  =  Cx  5Z  Ty7 
1=1  ,=1  V¥'  .=1  KV» 

which  implies  that 

‘mw  1 

c,  =  (Xi...  -  *„)  Ew  <7) 

,<=i  vy*. 

The  new  grid  intervals  along  the  x-axis  can  be  computed  from  (6)  as 


with  Cx  given  in  (7)  and  W{  given  in  (5),  and  the  coordinates  along  the  axis  are  given 
as 

t 

Xi  =  J2  Ax* 

k-l 

A  new  distribution  of  grid  points  along  the  y-axis  can  also  be  computed  using  the 
following  analogous  formulas. 

j 

Vi  =  5Z  &yk 

k=  1 
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r  -1  -1 

Jmax  ^ 

Cy  =  (yjmix  ~~  Vo)  ^  777~ 

J= 1  3 . 

Wj  =  1  +(u;~l)^Ml 

^max 

5  COMPUTATIONAL  RESULTS 

The  adaptive  rezoning  procedure  was  used  on  three  different  types  of  problems  in¬ 
volving  underwater  explosions.  The  procedure  was  most  successful  in  improving  the 
qualitative  rather  than  the  quantitative  nature  of  the  numerical  solution.  It  greatly 
reduced  the  oscillations  (or  ringing)  in  the  solution  without  the  addition  of  artificial 
viscosity.  However,  the  grids  were  still  not  fine  enough  in  the  neighborhood  of  shocks 
to  give  reliable  estimates  of  peak  pressures  and  impulse  values.  The  adaptive  grid 
algorithm  was  implemented  with  u  =  5  to  give  a  ratio  of  five  for  the  maximum- 
to-minimum  grid  spacing.  This  had  the  effect  of  reducing  the  grid  spacing  in  the 
neighborhood  of  the  shock  by  a  factor  of  three. 

The  gridding  scheme  was  capable  of  producing  extremely  fine  grids  near  the  shock, 
but  the  HULL  code  documentation  recommended  that  the  grid  aspect  ratio  not  exceed 
three.  This  ratio  was  exceeded  on  some  of  our  computations  with  no  noticeable  loss 
of  accuracy,  but  in  some  cases  erroneous  results  weie  obtained  when  the  aspect  ratio 
was  extremely  large.  The  weight  function  was  computed  with  cq  =  C\  =  1  and  C2  =  0, 
since  it  is  doubtful  that  the  second  derivatives  can  be  reliably  approximated  for  these 
types  of  problems.  Only  one  smoothing  iteration  was  used.  No  artificial  viscosity 
was  used  in  the  first  two  examples.  There  is  essentially  no  increase  in  core  storage 
with  the  adaptive  grid  since  the  code  is  already  written  with  variable  spacing  in  each 
coordinate  direction. 

The  first  example  is  the  computation  of  the  underwater  explosion  of  a  cylindrical 
charge  of  TNT.  This  was  solved  as  a  one-  dimensional  problem  with  a  500-by-'2  grid. 
The  radius  of  the  charge  was  7  cm  and  the  grid  extended  2500  cm  from  the  axis  of 
the  charge.  Figures  1  and  2  may  be  used  to  compare  the  solutions  at  t  =  0.1  seconds, 
computed  using  a  uniform  and  an  adaptive  grid.  The  location  of  every  tenth  cell  is 
indicated  along  the  top  borders  of  the  plots.  The  adaptive  grid  clearly  generated  a 
smoother  and  more  reclistic  pressure  profile;  however,  neither  peak  pressure  was  near 
the  theoretical  value  given  in  Cole[4].  In  this  particular  example,  the  peak  for  the 
uniform  grid  was  actually  closer  to  the  reference  value,  but  in  other  computations  the 
opposite  was  true.  The  adaptive  rezoning  increased  the  CPU  time  by  a  factor  of  two. 
However,  it  would  have  taken  three  times  as  many  grid  points  with  a  uniform  grid  to 
achieve  the  same  resolution  near  the  shock  wave. 

The  second  example  is  the  solution  of  an  axisymmetric  problem.  A  250-gram 
charge  of  TNT,  in  the  shape  of  a  cylinder  with  a  height  of  10  cm  and  a  radius  of  7 
cm,  is  detonated  underwater.  Two  plots  show  the  solution  computed  on  a  uniform 
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grid.  Figure  3  is  a  plot  at  an  early  time  when  the  shock  wave  is  near  the  charge,  and 
Figure  4  is  a  plot  at  a  later  time  when  the  shock  wave  has  reached  the  outer  boundary 
of  the  computational  region. 

Oscillations  in  the  solution  can  be  observed,  especially  near  the  axis  of  symmetry. 
The  solution  also  fails  to  develop  into  a  spherically  symmetric  solution  as  discussed 
by  Cole  [4].  Two  plots  are  also  included  for  the  solution  of  the  same  problem  on  an 
adaptive  grid.  Figure  5  displays  the  high  concentration  of  grid  points  near  the  charge 
that  is  needed  at  the  early  stage  of  the  computations  in  order  to  keep  oscillations  from 
initiating  and  propagating.  At  a  later  time,  as  indicated  in  Figure  6,  the  distribution 
of  grid  points  becomes  more  uniform  as  the  shock  front  moves  out  over  a  larger  section 
of  the  computational  region.  The  pressure  contours  and  velocity  profiles  are  also  much 
more  spherical  in  shape. 

The  final  example  demonstrates  the  application  of  adaptive  rezoning  to  the  solu¬ 
tion  of  problems  involving  several  materials.  The  explosive  is  a  spherical  charge  of 
TNT  weighing  66  kilotons.  The  TNT  is  lying  on  a  concrete  bottom  which  is  covered 
by  400  ft  of  water.  An  axisymmetric  solution  is  computed  with  a  column  of  air  ex¬ 
tending  to  a  height  of  200  ft  above  the  water  and  the  concrete  extending  100  ft  below 
the  water.  The  computational  region  is  truncated  at  a  distance  of  700  ft  from  the  axis 
of  symmetry.  Since  this  problem  took  a  lot  of  computer  time,  it  was  only  solved  with 
the  adaptive  rezoning.  Even  with  the  adaptive  grid,  no  reasonable  solution  could  be 
obtained  without  artificial  viscosity.  Figure  7  illustrates  the  solution  at  an  early  time. 
The  concentration  of  grid  points  at  the  axis  of  symmetry  was  especially  important  at 
this  stage  of  the  solution. 

Large  negative  pressures  were  produced  along  the  intersection  of  the  axis  of  sym¬ 
metry  and  the  concrete.  A  later  solution  is  plotted  in  Figure  8.  The  shock  wave  and 
the  fine  grid  region  have  moved  away  from  the  axis  of  symmetry.  The  full  horizon¬ 
tal  extent  of  the  grid  has  not  been  included  in  these  plots  so  that  the  details  of  the 
solution  can  be  more  easily  seen. 


6  CONCLUSIONS  AND  RECOMMENDATIONS 


It  has  been  demonstrated  that  adaptive  rezoning  can  significantly  improve  the  qual¬ 
ity  of  numerical  solutions  computed  by  the  HULL  code.  For  problems  involving 
explosions,  the  adaptive  grid  algorithm  can  be  automated  so  that  the  grid  points  are 
concentrated  in  regions  of  high  pressure  gradient.  As  with  all  short-term  projects,  we 
did  not  have  time  to  complete  all  that  we  would  have  liked.  The  code  still  does  not 
have  a  three-  dimensional  rezoning  capability.  The  rezoning  procedure  appears  to  be 
straightforward,  and  there  does  not  seem  to  be  any  problem  in  extending  it  to  three 
dimensions.  We  also  did  not  modify  the  code  so  that  several  rezoning  options  could 
be  used  simultaneously.  With  the  new  adaptive  method,  there  are  presently  nine  re¬ 
zone  options  in  the  HULL  code.  For  certain  types  of  problems  it  may  be  desireable  to 
use  more  than  one  type  of  rezoning.  For  example,  one  may  want  to  both  adaptively 
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redistribute  the  grid  points  and  expand  or  translate  the  physical  region  at  the  same 
time.  This  cannot  be  done  as  the  code  is  now  structured;  however,  the  rezone  option 
can  be  changed  between  successive  runs. 

The  construction  of  solution  adaptive  grids  is  limited  by  the  current  grid  restric¬ 
tions  in  the  HULL  code.  A  long-term  effort  should  be  made  to  lift  the  restriction  of 
having  to  compute  on  a  single  rectangular  grid  formed  by  the  Cartesian  product  of 
two  one-dimensional  grids.  Two  projects  could  be  undertaken,  either  of  which  would 
greatly  expand  upon  the  utility  of  the  current  code.  First,  the  differential  equations  of 
continuum  mechanics  can  be  formulated  in  terms  of  arbitrary  curvilinear  coordinates 
and  solved  on  a  curvilinear  coordinate  system.  This  would  require  a  reformulation 
of  the  finite  difference  algorithm  and  added  storage  for  additional  terms  and  variable 
coefficients  in  the  difference  equations.  However,  it  would  allow  the  user  to  fit  the 
boundary  lines  or  surfaces  to  the  true  charge  shapes  and  other  physical  boundaries 
in  the  problem.  The  grid  lines  could  also  be  placed  so  that  they  more  nearly  follow 
shock  waves  in  the  solution.  Of  course,  a  true  two-dimensional  rezoning  could  be 
done  on  this  type  of  grid  structure. 

A  second  enhancement  of  the  code  would  be  the  ability  to  compute  solutions  on 
multiple  arrays  of  grid  points.  Each  array  of  points  would  correspond  to  a  rectan¬ 
gular  block  and  each  block  would  have  its  own  grid  spacing  in  each  direction.  In  an 
adaptive  implementation,  rectangular  sub-blocks  of  a  global  grid  could  be  cut  out  and 
replaced  by  blocks  with  more  grid  points.  When  using  rectangular  coordinates,  much 
of  the  HULL  code  would  remain  unchanged.  A  new  main  program  would  be  needed 
which  would  organize  the  sweeps  over  each  rectangular  block  of  points  and  transfer 
information  between  blocks  at  each  time  step.  In  the  case  of  unequal  grid  spacing 
at  block  interfaces,  interpolation  schemes  would  have  to  be  developed  which  preserve 
the  conservation  properties  of  the  numerical  solution.  The  initial  effort  could  begin 
with  either  of  these  projects.  If  both  are  eventually  completed,  the  result  would  be  a 
code  which  could  treat  multiple  blocks  of  curvilinear  grid  systems  of  any  size  or  shape 
and  any  grid  point  density.  The  basic  code  would  still  compute  on  rectangular  arrays 
so  that  there  would  be  much  greater  freedom  in  positioning  grid  points  without  the 
overhead  costs  associated  with  completely  unstructured  grids. 
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8  APPENDIX  A 


This  appendix  will  explain  how  to  use  the  adaptive  rezone  option  and  the  grid  plotting 
capability  which  have  been  added  to  the  HULL  code.  A  rezone  option  is  called  by 
specifying  a  value  for  the  parameter  REZONE  in  the  HULL  Euler  module  main 
routine.  For  the  new  adaptive  option,  set 

REZONE  =  9 

after  the  keyword  INPUT  on  the  input  file.  All  parameters  in  the  rezone  algorithm 
have  been  set  to  values  which  appear  to  give  the  optimal  grid  for  the  problems  which 
have  been  considered  in  this  report.  They  are  set  in  the  following  data  statement, 
which  can  be  changed  using  the  SAIL  edit  facilities. 

DATA  OMEGA ,  COEFO,  COEF\,COEF2,  NOSMOO/5.,  1.,  1.,  0.,  1/ 
Relating  these  variables  to  the  notation  in  Section  2  of  this  report, 


OMEGA 

=  U ) 

COEFO 

=  Co 

COEFl 

=  Cj 

COEF2 

=  c2, 

and  NOS  MOO  is  the  number  of  smoothing  iterations  applied  to  the  function  P. 

Plots  of  the  computational  grid  can  be  made  by  running  program  PULL.  The 
code  for  producing  the  grid  plots  was  inserted  with  only  a  few  changes  in  the  PULL 
program.  A  plot  of  the  grid  is  the  default  option  when  making  any  contour  plot. 
Thus,  pressure  contours  and  the  grid  would  be  plotted  by  inserting  the  keyword 

PCONT 

in  the  data  file.  The  parameter  used  to  suppress  the  printing  of  cell  indices  is  also 
used  to  suppress  grid  plotting.  For  example,  to  plot  pressure  contours  without  a  grid, 
use 

PCONT  (  NOCELLS  ). 

Since  the  grid  is  generated  during  a  contour  plot,  it  is  not  possible  to  generate  a 
separate  plot  of  only  the  grid.  It  is  also  not  possible  to  print  cell  indices  on  the 
border  of  a  contour  plot  without  drawing  the  grid. 
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Radius  m 

Time  10.001  ms  Cycle  2150  Problem  1001.0000 

UNDERWATER  EXPLOSION 


Figure  1.  Pressure  field  from  the  detonation  of  an  infinite  cylinder  of  TNT 
explosive,  t  =  10  msec,  computed  on  a.  uniform  grid. 
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dynes/cm*  (xIO7) 


Figure  2.  Pressure  field  from  the  detonation  of  an  infinite  cylinder  of  TNT 
explosive,  t  =  10  msec,  computed  on  an  adaptive  grid. 
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Figure  3.  Calculated  values  for  a  finite  cylindrical  charge  weighing  250  grams,  t  = 

200  n  sec. 
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Figure  4.  Calculated  values  for  a  finite  cylindrical  charge  weighing  250  grams, 

t  =  1  msec. 
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Figure  5.  Calculated  values  on  an  adaptive  grid,  t  =  200  /x  sec. 
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Figure  6.  Calculated  values  on  an  adaptive  grid,  t  =  1  msec. 
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Figure  7.  Explosion  of  66  kilotons  of  TNT  in  400  ft  of  water,  t  =  20  msec. 
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Figure  8.  Explosion  of  66  kilotons  of  TNT  in  400  ft  of  water,  t  =  100  msec. 
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